      subroutine createbp(ntie,ipkon,kon,
     &  lakon,nslav,nmast,imast,nfslav,
     &  islavsurf,imastsurf,pmastsurf,itiefac,contribution)
!
!     compute the Bp matrix for contact problems
!     Author: Li, Yang
!
      implicit none
!
      character*8 lakon(*)
!
      logical contained
!
      integer ntie,ipkon(*),kon(*),iflag,i,j,k,l,
     &  ifaceq(8,6),ifacet(6,4),ifacew(8,5),
     &  indexe,nope,islavsurf(2,*),
     &  imastsurf(*),itiefac(2,*),ifaces,nelems,jfaces,ifacem,nelemm,
     &  jfacem,mint2d,indexf,m,nopes,nslav,nmast,imast,
     &  nfslav,iq,ii,ig,nodem
! nslav: num of slave nodes
! nmast: num of master nodes
! imast: all master nodes
! nfslav: num of slave faces
!
      real*8 pmastsurf(2,*),
     &  ets,xis,etm,xim,weights,xl2(0:3,8),xsj2(3),
     &  shp2(4,8),xs2(3,2),xld(0:3,8),xsjd(3),shpd(4,8),xsd(3,2),
     &  contribution
!
      include "gauss.f"
!
      data ifaceq /4,3,2,1,11,10,9,12,
     &            5,6,7,8,13,14,15,16,
     &            1,2,6,5,9,18,13,17,
     &            2,3,7,6,10,19,14,18,
     &            3,4,8,7,11,20,15,19,
     &            4,1,5,8,12,17,16,20/
      data ifacet /1,3,2,7,6,5,
     &             1,2,4,5,9,8,
     &             2,3,4,6,10,9,
     &             1,4,3,8,10,7/
      data ifacew /1,3,2,9,8,7,0,0,
     &             4,5,6,10,11,12,0,0,
     &             1,2,5,4,7,14,10,13,
     &             2,3,6,5,8,15,11,14,
     &             4,6,3,1,12,15,9,13/
!
      do j = 1,nslav
         do k = 1,nmast
            contribution = 0.d0
            nodem = imast(k)
            do i = 1,ntie
            do l = itiefac(1,i), itiefac(2,i)
               ifaces = islavsurf(1,l)
               nelems = int(ifaces/10)
               jfaces = ifaces - nelems*10
!
! Decide the max integration points number, junt consider 2D situation 
! 
                if(lakon(nelems)(4:5).eq.'8R') then
                   mint2d=1
                elseif((lakon(nelems)(4:4).eq.'8').or.
     &             (lakon(nelems)(4:6).eq.'20R')) then
                   mint2d=4
                elseif(lakon(nelems)(4:4).eq.'2') then
                   mint2d=9
                elseif(lakon(nelems)(4:5).eq.'10') then
                   mint2d=3
                elseif(lakon(nelems)(4:4).eq.'4') then
                   mint2d=1
                endif
!
!         treatment of wedge faces
!
                if(lakon(nelems)(4:4).eq.'6') then
                   mint2d=1
                   if(jfaces.le.2) then
                      nopes=3
                   else
                      nopes=4
                   endif
                endif
                if(lakon(nelems)(4:5).eq.'15') then
                   if(jfaces.le.2) then
                      mint2d=3
                      nopes=6
                   else
                      mint2d=4
                      nopes=8
                   endif
                endif
!
                indexf = islavsurf(2,l)
                do m = 1,mint2d
                   ifacem = imastsurf(indexf + m)
                   nelemm = int(ifacem/10)
                   jfacem = ifacem - nelemm*10
!
                   indexe = ipkon(nelemm)
!
                   if(lakon(nelems)(4:4).eq.'2') then
                      nope=20
                      nopes=8
                   elseif(lakon(nelems)(4:4).eq.'8') then
                      nope=8
                      nopes=4
                   elseif(lakon(nelems)(4:5).eq.'10') then
                      nope=10
                      nopes=6
                   elseif(lakon(nelems)(4:4).eq.'4') then
                      nope=4
                      nopes=3
                   elseif(lakon(nelems)(4:5).eq.'15') then
                      nope=15
                   elseif(lakon(nelems)(4:4).eq.'6') then
                      nope=6
                   endif
!     
!     treatment of wedge faces
!     
                   if(lakon(nelems)(4:4).eq.'6') then
                     if(ig.le.2) then
                        nopes=3
                     else
                        nopes=4
                     endif
                   endif
                   if(lakon(nelems)(4:5).eq.'15') then
                     if(ig.le.2) then
                        nopes=6
                     else
                        nopes=8
                     endif
                   endif
!     
                   contained = .FALSE.
                   if((nope.eq.20).or.(nope.eq.8)) then
                     do ii=1,nopes
                        iq=kon(indexe+ifaceq(ii,jfacem))
                        if(iq.eq.nelemm) then
                            contained=.TRUE.
                            exit
                        endif
                     enddo
                   elseif((nope.eq.10).or.(nope.eq.4)) then
                     do ii=1,nopes
                        iq=kon(indexe+ifacet(ii,jfacem))
                        if(iq.eq.nelemm)  then
                            contained=.TRUE.
                            exit
                        endif
                     enddo
                   else
                     do ii=1,nopes
                        iq=kon(indexe+ifacew(ii,jfacem))
                        if(iq.eq.nelemm) then
                            contained=.TRUE.
                            exit
                        endif
                     enddo
                   endif
!     
                   if(.not.contained)  cycle
!
                   if((lakon(nelems)(4:5).eq.'8R').or.
     &               ((lakon(nelems)(4:4).eq.'6').and.
     &               (nopes.eq.4))) then
                      xis=gauss2d1(1,m)
                      ets=gauss2d1(2,m)
                      weights=weight2d1(m)
                   elseif((lakon(nelems)(4:4).eq.'8').or.
     &                   (lakon(nelems)(4:6).eq.'20R').or.
     &                   ((lakon(nelems)(4:5).eq.'15').and.
     &                   (nopes.eq.8))) then
                      xis=gauss2d2(1,m)
                      ets=gauss2d2(2,m)
                      weights=weight2d2(m)
                   elseif(lakon(nelems)(4:4).eq.'2') then
                      xis=gauss2d3(1,m)
                      ets=gauss2d3(2,m)
                      weights=weight2d3(m)
                   elseif((lakon(nelems)(4:5).eq.'10').or.
     &                   ((lakon(nelems)(4:5).eq.'15').and.
     &                   (nopes.eq.6))) then
                      xis=gauss2d5(1,m)
                      ets=gauss2d5(2,m)
                      weights=weight2d5(m)
                   elseif((lakon(nelems)(4:4).eq.'4').or.
     &                   ((lakon(nelems)(4:4).eq.'6').and.
     &                   (nopes.eq.3))) then
                      xis=gauss2d4(1,m)
                      ets=gauss2d4(2,m)
                      weights=weight2d4(m)
                   endif
!
                   iflag = 2
                   if(nopes.eq.8) then
                     call shape8q(xis,ets,xl2,xsj2,xs2,shp2,iflag)
                   elseif(nopes.eq.4) then
                     call shape4q(xis,ets,xl2,xsj2,xs2,shp2,iflag)
                   elseif(nopes.eq.6) then
                     call shape6tri(xis,ets,xl2,xsj2,xs2,shp2,iflag)
                   else
                     call shape3tri(xis,ets,xl2,xsj2,xs2,shp2,iflag)
                   endif
!
                   xim = pmastsurf(1,indexf+m)
                   etm = pmastsurf(2,indexf+m)
!
                   if(nopes.eq.8) then
                     call dualshape8q(xim,etm,xld,xsjd,xsd,shpd,iflag)
                   elseif(nopes.eq.4) then
                     call dualshape4q(xim,etm,xld,xsjd,xsd,shpd,iflag)
                   elseif(nopes.eq.6) then
                     call dualshape6tri(xim,etm,xld,xsjd,xsd,shpd,iflag)
                   else
                     call dualshape3tri(xim,etm,xld,xsjd,xsd,shpd,iflag)
                   endif
                   contribution=contribution+shp2(4,m)*shpd(4,m)*
     &                   weights*dsqrt(xsj2(1)**2+xsj2(2)**2+xsj2(3)**2)
                enddo
            enddo
        enddo
      enddo
      enddo
!
      return 
      end